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■ We study analytically the coarse and fine-grained distribution function 



ABSTRACT 



established by the self-similar infall of collisionless matter. We find this func- 
tion explicitly for isotropic and spherically symmetric systems in terms of 
cosmological initial conditions. The coarse-grained function is structureless 

£3 ■ 

and steady but the familiar phase space sheet sub-structure is recovered in 

Cd ■ 

the fine-grained limit. By breaking the self-similarity of the halo infall we are 
able to argue for a central density flattening. In addition there will be an edge 
steepening. The best fitting analytic density function is likely to be provided 
by a high order polytrope fit smoothly to an outer power law of index —3 for 
isolated systems. There may be a transition to a —4 power law in the outer 
regions of tidally truncated systems. 

Since we find that the central flattening is progressive in time, dynami- 
cally young systems such as galaxy clusters may well possess the NFW type 
density profile, while primordial dwarf galaxies for example are expected to 
have cores. This progressive flattening is expected to end either in the non- 
singular isothermal sphere, or in the non-singular metastable polytropic cores; 
since the distribution functions associated with each of these arise naturally 
in the bulk halo during the infall. We suggest based on previous studies of 
the evolution of de-stabilized polytropes, that a collisionless system may pass 
through a family of polytropes of increasing order, finally approaching the 
limit of the non-singular isothermal sphere, if the 'violent' collective relaxation 
is frequently re-excited by 'merger' events. Thus cD galaxies and indeed all 
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bright galaxies that have grown in this fashion, should be in polytropic states. 
Our results suggest that no physics beyond that of wave-particle scattering is 
necessary to explain the nature of dark matter density profiles. However this 
may be assisted by the scattering of particles from the centre of the system 
by the infall of dwarf galaxies, galactic nuclei or black holes (e.g. Nakano and 
Makino, 1999), all of which would restart the pure dynamical relaxation. 

Key words: methods: analytical-gravitation-accretion- distribution func- 
tions. 



1 INTRODUCTION 

In a previous paper (Henriksen and Le Delliou, 2002 = paper I) a novel analytic method 
for studying the (mass) Distribution Function (DF) of a collisionless self-gravitating system 
of particles was introduced. The essential idea is to carry out an expansion of the coupled 
Boltzmann-Poisson system in terms of the phase space resolution. Normally the zeroth 
order in the expansion corresponds to the coarsest-grained view, while the higher order 
terms correspond to less crude coarse-graining. The variable resolution of the system phase 
space is achieved by means of appropriate non-canonical coordinate transformations, these 
being chosen to contain a parameter whose magnitude adjusts the minimum phase space 
volume resolved. A fully relaxed system can be expected to look the same at different levels 
of coarse-graining and so this method permits the construction of equilibrium distribution 
functions by determining the conditions under which all structure at any level of coarse- 
graining vanishes. A fine-grained expansion also exists wherein an expansion about the 
coarse-grained structureless DF in a series of increasing resolution recovers the phase-space 
structure known to exist in the infall models. 

It should be noted that the value 1 of the expansion parameter a corresponds to the 
phase space volume that is used to define a smooth DF. For example, were we applying the 
technique to a clumped dark matter system or to a system of globular clusters, there would 
be many individual objects in the phase space volume used in the definition. Larger volumes 
correspond to an even coarser graining and so terminating the series in 1/a at zeroth order 
is a means of confirming the smoothed nature of the original definition. We obtain fine 
structure by expanding in values of a < 1, as is shown below. This provides assurance that 
the resolution expansion can recover essential features of the DF. 

One knows however that the strict self-similarity in the infall must be broken near the 
centre of the system (paper I and below, where our series may be seen not to converge as 
the potential gradient diverges). This leads us to study the explicit time dependence of the 
system as the direct means of breaking the self-similarity (Carter and Henriksen, 1991) and 
this predicts a central flattening in time. In Merrall and Henriksen (2003) a de-stabilized 
polytrope was found to relax towards a gaussian DF with dispersion that increased in time. 
That is the 'temperature' of the system increased until a new equilibrium was achieved. This 
effect is also observed in our analytic study of the central relaxation. 

The technique for generating the structureless DF is best employed for steady state 
systems since one expects relaxed systems to be steady. However there is a particular case of 
time dependence that is of importance since it corresponds to a possible growth mode of dark 
matter halos. This is the well studied self-similar infall (e.g. Henriksen and Widrow, 1999 
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and references therein, Le Delliou and Henriksen, 2003), which is known to arise inevitably 
for radial orbits and can be assumed as a general model. In the appropriate variables this 
growing mode of the system can be treated as steady, and in fact we find that the DF's 
constructed in this way are absolutely steady (indeed they may be constructed by this 
double requirement; Henriksen and Widrow, 1999). Thus they represent equilibrium states 
of the system, at least until they are subsequently disturbed. 

A first study of these systems was conducted in paper I. In particular that paper applied 
the method to dark matter halos in spherical symmetry to show that the self-similar density 
profile had to flatten at the centre of the system. It was found that the flattening would be 
progressive, with the effective power law of the density profile passing through the NFW 
(Navarro, Frenk and White, 1996) value of —1 on the way to even flatter values. This 
conclusion seems now to be in agreement with the numerical simulations (e.g. Power et al., 
2003), which fact offers confirmation to both methods. 

The discussion of isotropically elliptical orbits in paper I is unfortunately flawed by the 
omission of a factor (i in the last term of equation (80) of that paper that was carried 
through to the end of that section. Since this term vanishes when j — 3, the discussion of 
that case remains valid. The general formulae are corrected in the appendix to this paper. 

In the bulk of the present paper we present a conclusive application of our ideas to an 
isotropic, spherically symmetric collisionless dynamical system. This case is likely to be a 
general ultimate equilibrium since theory, as elaborated e.g. by Lynden-Bell (1967), Aly 
(1989), Nakamura (2000), suggests that the ultimate relaxed state of a collisionless system 
should be spherically symmetric with an isotropic Gaussian DF in energy. However our 
techniques may be applied to a general system and this is discussed in a mathematical 
appendix. This allows a discussion of axial symmetry and of the dependence on phase-space 
dimensionality in general, but the physical conclusions are not different from those presented 
in the text when considerations of regularity, symmetry and entropy are applied. 

The mass distribution of the dominant dark matter in many dwarf elliptical galaxies 
is best fit by that of the non-singular isothermal sphere (e.g. Cote, Carignan & Freeman, 
2000; Begum, Chengalur and Hopp, 2003; Begum and Chengulur, 2004). On the other hand 
XMM and Chandra studies of hydrostatic hot gas in galaxy clusters (e.g. Buote, 2004; Pratt 
and Arnaud, 2003, 2002) show in most cases that the NFW profile gives a good fit to the 
total mass distribution down to about .03-. 05 r 2 oo or about 40- 70 kpc for the low mass 
cluster Abel 1983. However the fact that the virial radius r2oo is only about 4 times the core 
radius suggests that perhaps the 'real' core has not been detected. Based on the results of 
the present paper we would expect such a low mass cluster to have a core relatively large 
compared to the high mass clusters. 

Sand et al. (2002) and Gavazzi et al. (2003) have however studied the cluster MS 2137.3 
-2353 using gravitational lensing of background sources and conclude that the central power 
law is as flat as —0.35. This has been 'explained' in terms of infalling baryonic clumps (El- 
Zant et al. 2004) but we find that such flattening occurs naturally during the relaxation of 
the dark matter infall and that there is no need for such extrnal influences. It is undoubtedly 
the case that such 'clumpy' infall or merging does have a flattening influence (e.g. Nakano 
and Makino, 1999), and we note that this happens automatically in the self-similar infall 
since shells that fall in later have more mass (for j < 3/2). 

A subsequent paper will attempt to convert our predictions as a function of dynamical 
age to a mass dependence. In the present work we wish principally to justify the method 
and the formal results. In the next section we present our results in the context of an 
isotropic, spherically symmetric halo, since these correspond to our general results and are 
more accessible. The appendix will provide more mathematical detail for the general case. 
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2 ISOTROPIC, SPHERICALLY SYMMETRIC, HALOS 

We present in this section the formulae that are common to both the general case treated 
in the appendix and to the spherically-symmetric isotropic case treated here. 



2.1 Mathematical formulation 

As usual we employ a mass distribution function and here it is taken to be in six phase space 
dimensions so that / = /(r, v; t). Thus we have to solve 

$/ + *v/ + -.v„/ = o, (1) 

and 

V 2 $ = AnGp, (2) 
where 

P = J f d 3 v. (3) 

For some applications it is possible to add an external potential to the 'internal' $. This 
could be due to a central point mass or to a distribution of baryons. 
In spherical polar coordinates we must remember that 

dv = vJ + 4 d ^ V\ COt(6>) - VrVQ - d e $ VjjVr + vo cot(fl)) + 4gy^g 

dt r ' r r 

The variable transformations that render the problem stationary during the self-similar 
infall are: 

R = re~ ST , Y = ve~^ T , * = <$>e 2 ^ T , 
P = /e (3f-i)aT Q= pe (2a)r ? e «r = af . (5) 

The angles # and 0, being dimensionless, do not transform. 
The equations in these transformed variables become 

WtP _ «. 1)P+( L^ + (l^.(!. 1)y| .M) v 

a a a \ ai? a a J 

\ an a J 

an a J 

and 

= AtvG 0, (7) 

where 

= 1 Pd 3 Y", (8) 

and the R indicates the scaled radial coordinate (the angular coordinates are unchanged). 
For the isotropic and spherically symmetric limit that we study in this paper, we set 

P = P(Y 2 ,R;t) # = W(i2;t), (9) 

in the previous equations, where 
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Y* = Y* + Yl = Y* + Y? + Y*. (10) 
Then from a direct substitution into equation (6) we obtain 

tdtP - (3- - l)P + - -R)d R P -(-- l)Yd Y P - (ld R *)&)dyP = 0, (11) 

a a a a a Y 

and of course equation (7) becomes simply 

The coarse-graining or fine-graining is carried out in terms of the parameter a since 
it is easy to show from the transformations of equation (5) that the physical phase space 
element is proportional to the factor exp (6- — 3)aT multiplying the phase space element 
in transformed variables. Thus for a fixed self-similar index (^) we may increase or decrease 
the physical phase space volume (relative to the reference value at a = 1) corresponding to a 
given set of transformed variables by increasing or decreasing a respectively. This is true for 
j < 2. This limiting value 2 renders the transformation canonical so that no transformation 
of the phase space volume is effected. We can nevertheless approach this limit arbitrarily 
closely; but in any case, as discussed in the appendix, an expansion in powers oil/ a amounts 
to an expansion in time resolution and therefore the lowest order effectively chooses the 
steady state DF . 

The coarse graining series is taken as 

P = P Q + P l a- 1 + p 2 a~ 2 + ... (13) 
while the fine graining series is assumed in the form 

P = P + P- 1 a + P^ 2 a 2 + ... (14) 

We apply the coarse-graining expansion in the next sub-section and follow it in a sub- 
sequent sub-section with a discussion of the results. A final sub-section in this part applies 
the fine-graining expansion to show that the 'sheet' phase-space structure may be recovered 
in this way. 



2.2 Coarse Graining 

We substitute the coarse graining expansion into equation (ll)and obtain at zeroth order 
the equation 

A A A 

-td t P + (3- - 1)P + (-)Rd R P + (- - l)Ydy P = 0, (15) 
a a a 

while the second equation at zeroth order is simply, from equation (12), 

Using the method of characteristics in these latter two equations we obtain the solution 

Po = Poo(()e-^- 1)s , O = IooR-^\ t = t o (0e~ s , 

^ o = _ 7 i?2(i-f) )C =^-1)^ s =^\nR, (17) 

o 

where 

7 = (47rG/ 00 )/(2(3-2^)(|-l)), (18) 
and 
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loo = J Poo(C) d 3 0- (19) 

In the preceding formulae, we have assumed that all particles cross R — 1 at s — 
(although not necessarily at the same time of course) with the scaled velocity Y equal to 
(. In general we may identify s with time through aT = ln(o;t ) — s. However there is 
no assumption of self-similarity until the explicit dependence on T or t is suppressed. We 
continue for the moment by assuming self-similarity so that d t = 0, but we will recall below 
that it is an explicit means of breaking the self-similarity. 

In the special case having f = 1 the above formulae apply except that 

^ = 47rG/ 00 ln(i?). (20) 

It should also be noted that in that case Y = v = (. 

The first order coarse-grained Vlasov equation now becomes 

5 5 5 Y 

(3- - + -Rd R P 1 + (- - l)Yd Y P l - Y R d R P Q + (-^)d R V dyP = 0, (21) 
a a a Y 

from which we see that the characteristics do not change with coarse grained order and in 
addition, after a short calculation which takes care to differentiate before evaluating on the 
characteristics, that 

+ <s£ - l)Pi = e-^C- - 1)0, (« - 2 7 /0^f - pP.) • (22) 

It is readily seen that all higher order terms in the expansion will vanish if the first order 
correction vanishes (apart from the homogeneous part which has the same form as P OQ and 
can be absorbed at this order). Thus we obtain the condition 

(C-2 7 /C)^-|^|p oo = 0, (23) 

to obtain an exact smoothed DF independent of the degree of coarse-graining. 
When f = 1 the first order equation becomes 

^ +2 P 1 = _^( 2Poo + !£^!?M (24) 



ds R 3 \ v dv 

and so the perfectly coarse-grained condition takes the form 

^ + TTT P " = °- < 25 > 

d v 47rGi OD 
Equations (23) and (24) integrate to give respectively 

Poo = C\e 00 \^f\ (26) 
and 

P OQ = C e '^/^Gi 00 )_ (27) 
Here we have set 

.a, 3 — % _ N 

*i> - Kf^ry < 28 > 

and 

e 00 = C 2 /2 - 7- (29) 
This latter quantity is related to the energy E = $ + v 2 /2 by 

E = e 2 ^ 5T R 2 ^e 00 . (30) 
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Equation (26),together with equations (17) and (5), yields an exact, steady (in agreement 
with Jeans' theorem), coarse-grained DF in the polytropic form (E < for convergence) 

f = C\E\ p W. (31) 
A similar procedure in the case f = 1 yields the Gaussian 

/ = Ce~^, (32) 

where£ = v 2 /2 + inGI 00 ln(r). 

We conclude that in this phase space geometry, self-similar infall combined with the 
associated collective relaxation leads to a smoothed DF in the above form. It is generally a 
polytrope of order 

n=j/(j-l), (33) 

since the Gaussian form is a limiting polytrope of infinite order which is attained as j — > 1+. 
We expect only the cases f > 1 since otherwise the polytropes are linearly unstable according 
to the Doremus-Feix-Baumann theorem (DFBT; e.g. Binney and Tremaine, 1987). 



2.3 Discussion of Coarse Grained Results 

Although our unique results (31) and (32) are found here for a spherically symmetric isotropic 
system, we emphasize that similar results may be found using this method for other sim- 
ple systems such as the ID motion of collisionless planes, purely radial orbits in spherical 
symmetry and anisotropic spherical symmetry. They are similar in the sense that a similar 
sequence of gaussian polytropes is found, but they differ in index depending on the assumed 
phase space dimensions. This dependence is in accord with the numerical conclusions of 
Moutarde et al. (1995), Teyssier et al. (1997) and Torman et al. (1997). In addition we 
show in the appendix that this same sequence is accessible starting from the most general 
geometry in phase space. 

We should note that Merrall and Henriksen (2003) have shown by several independent 
numerical simulations that even linearly stable polytropes are in fact metastable, and that 
they collapse after sufficient cooling towards a (lowered) Gaussian distribution, whose tem- 
perature increases in time over the initial free-fall value. Given this result, and given that 
theory as elaborated e.g. by Lynden-Bell (1967), Aly (1989), and Nakamura (2000) suggests 
that the ultimate DF of a relaxed collisionless system should be a spherically symmetric, 
isotropic Gaussian, at least sufficiently far from all boundaries; we suggest that the collec- 
tive relaxation proceeds through a polytropic sequence of distribution functions as above. 
Moreover since Merrall and Henriksen (ibid) also find that the relaxation 'stalls' if there is 
too little initial symmetry, we suggest further that a system may spend considerable time in 
a metastable polytropic state before being re-excited by some external encounter and thus 
being provoked to continue its evolution towards a Gaussian (energy dissipation) or towards 
a lower order polytrope if the energy exchange is positive. 

We should also recall in this connection that 'initially steep' (that is % > 1 in the 
self-similar regime) halos were shown by Henriksen and Widrow (1997) to be isolated from 
their surroundings by an external saddle point in the family of possible orbits, and this 
was suggested as the origin of the 'memory' of initial conditions. There is no such point for 
f < 1 and these evolve rapidly towards f = 1- However although suggestive, these results 
are strictly only for radial orbits. 

In the present case the characteristics of equation (11) do not in general admit the same 
analytic treatment. But in the limiting case of j = 2 this treatment is in fact possible, and 
one finds, following the procedure of Henriksen and Widrow (1997), that (setting a — 1 for 
convenience in this argument) the equivalent one dimensional problem in R becomes 
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d*R _ _cM Ml 

ds 2 dR ' 1 ' 

where * e // = * - -R 2 / 8 - J 2 /{2R 2 ) and Fj_ = J/i? with J constant. Thus for negative 
energy particles there is always an outer turning point just as for purely radial orbits, and 
in addition there is the inner turning point required by non-zero angular momentum. Thus 
the isolation by critical points is also present in this limit, which is of course initially steep. 
Although we can not establish this behaviour for all steep indices this analytic argument adds 
weight to the numerical evidence for meta-stability cited above. Such behaviour suggests 
that the polytropic index will be preserved during the single halo infall and that subsequent 
disturbances are required to break the meta-stability. 

Since no finite system can realize either a global polytrope of n > 5 or a Gaussian sphere, 
we expect that the outer parts will be dominated by the edge effect power law r~ 3 (perhaps 
r~ 4 for tidally truncated cases). In the central bounded regions higher order polytropes are 
more likely. 

But it is the density profile which we must address next. This follows from equations 
(17) and (5) unequivocally as 

P = IooT^. (35) 

Moreover the similarity class may be expressed in terms of the primordial density per- 
turbation from which the halo is born (Fillmore and Goldreich, 1985; Bertschinger, 1985; 
Henriksen and Widrow, 1999) if this is assumed to have been a power law of the form 
5p oc r~ e . One finds 

^ = 3e/2(e+l), (36) 
and so the initial polytropic DF might be expected to have by equation (33) the order 

» = ^, (37) 
for e > 2. 

Thus formally for e > 2 we deduce a power law density profile in the bulk of the halo that 
is approximately p oc r -2 . This implies a singular density expression of both the polytropic 
and the gaussian distribution functions. However it can only apply where strict self-similarity 
is reasonable, which is always in an intermediate region far from boundaries. The polytropic 
and gaussian distribution functions are compatible with non-singular density distributions as 
well as with these singular, mathematically isolated, forms and so we must look for evidence 
that the system will tend to evolve toward the non-singular family. Henriksen and Widrow 
(1999) have shown that near the edge of such a system one can expect the 'Keplerian' power 
law p oc r -3 . They were able to give a self-consistent analysis for the region outside a fixed 
mass, but the simple version of their argument consists of taking e large as might be expected 
near the edge of the system, which gives j = 3/2 and hence the Keplerian power law. 

So we already expect the inverse square law, to turn over at the edge to become an 
inverse cube law. This accords qualitatively with the simulations (e.g. Navarro, Frenk and 
White, 1996) except at the centre of the system. But it is certainly at the centre of the 
system where we must break the coarse-graining self-similarity since the expansion will not 
converge where the potential gradient becomes infinite. We can only hope to approach this 
region asymptotically by breaking the self-similarity in an effort to find a tendency to flatten. 
This we do in a subsequent section although the essential idea was already discussed in paper 
I. 

Our tentative conclusion is that we should select from the family of polytropes, at least 
in the central regions to find a continuous prediction for the density profile. We believe it to 
be a high order polytrope that is fit smoothly to an outer edge effect r~ 3 law. For example 
as is shown in figure (1) a polytrope of order 8, for which j — 8/7 will connect smoothly a 



Coarse- graining the distribution function of cold dark matter II 9 



r 



-1 -0.5 0.5 1 




Figure 1. The left panel shows the log-log density profile of an n = 8 polytrope together with an r~ 3 power law (dotted). The 
constants in the power law have been chosen so that it fits continuously to the polytrope at r = 7.706. This is just beyond the 
point where the mass equals f and is the point where the polytropic slope is 3. The right panel shows the slope of the polytrope 
as a function of log r. 



region of slope equal to —3 to a central core. As might be expected, even the non-singular 
isothermal sphere(see e.g. Binney and Tremaine,ibid, p230 ) will join a region of slope nearly 
equal to —3 to a central core. The polytropic region must terminate in any case to keep the 
mass finite, when n > 5, and this edge effect is reponsible for the r~ 3 law. Strictly the outer 
slope should be slightly steeper than this in order to render the mass finite. Our suggested 
fit to the density profiles of dark matter halos replaces the phenomenological parameters of 
the NFW fit by the index of the central polytrope. Since this in turn may depend on the 
age of the system, we expect universal fits only for the oldest systems. 

One should note however that there are other physical effects that may fix the similarity 
class in addition to the mass distribution, such as a fixed specific angular momentum, for 
which j attains the limiting value of 2. This is because together with Newton's constant G 
such a quantity determines a density profile 

j 2 i 

P = V-q~v ( 38 ) 

where fj, is a positive numerical constant. We see from equation (36) that since we expect e to 
be > 2, there is no value that gives f = 2. Thus we must abandon setting the initial density 
profile independently of a limiting specific angular momentum J. A halo dominated by a 
constant specific angular momentum could be relevant to halos formed from material that 
has been spun up tidally (presumably that which has fallen in from larger spatial scales), 
and indeed tidally truncated. The outer regions of globular clusters and of elliptical galaxies 
may reflect this behaviour in their density profiles, and in fact any dark matter halo must 
ultimately become steeper than r~ 3 , probably due to tidal effects. Our argument requires 
the region of tidal truncation to have a density profile p oc r~ 4 . 

Because our transformation to stationary coordinates is canonical in this case, no coarse 
graining is actually achieved as a — > oo. However this should not affect the result that is 
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stopped at zeroth order which yields the steady polytrope, since as discussed in the appendix 
the expansion in 1/a can also be understood as a progressive time averaging. 

Finally in this section we note that we can determine the positive constant C in equation 
(32) by introducing spatial and density scales, and thus completely determine that result in 
physical terms. We introduce a fiducial density and a fiducial radius p Q , r Q , so that from the 
Poisson equation we have 4nGp = k/r 2 and hence 

k = AuG Po r 2 . (39) 

Then from the integral of the DF in (32) over velocity space (making the usual extention to 
infinite velocities) we have 

c = <!> 3 ' 2 x skf x 7T < 40 > 

In the next sub-section we turn attention to the fine graining limit in order to demonstrate 
that we can gradually recover the known phase space structure (e.g. Henriksen and Widrow, 
1999; Merrall and Henriksen, 2003). 



2.4 Fine Graining 

In this section we proceed by using the same transformation as in the previous section, but 
now the expansion is of the form 

P = P + aP_i + a 2 P_ 2 + . . . , (41) 

where P Q is the steady coarse-grained function found in the previous section, see equations 
(31,32), and the quantity \1/ G is the corresponding potential. Recalling the self-similar form 
(no explicit time dependence) of equation (15), this expansion yields at zeroth order in a 

W _ 1 _ ( ^ )( ^ )a ,P. 1 = ( ^L i)( ^ )a ,P . (42) 

From the characteristics of this equation we have 
dR dY Y R dm 

77 " Yr > 77 ~ ~ ( T K 71f j ' 
dP-i , dj^ ( R r 2 

— = ( 7nr )( 7 )arPo ' T + * o = e - (43) 

The scaled energy e is constant on the fine-grained characteristic and it takes the general 
form 

e = e 00 R^~f\ (44) 

with m o either as in equations (17, 18) for j > 1 or as in equation (20) for f = 1- The coarse 
grained zeroth order must also be adjusted accordingly when j — 1, becoming P OQ from 
equation (27) times R~ 2 . Otherwise f = 1 may be substituted directly into the equations 
above. The zeroth order coarse-grained characteristic ( is also as in equation (17). 

Proceeding now with the general case of equation(31), we calculate from the appropriate 
member of equations (43) that 

We must now use Poisson's equation in the form 

^•*7n(* d -iir)' (46) 
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with 

©-! = / P-i d 3 Y = R^ 1 -^ J P^d 3 (, (47) 

to write an equation for after which we find 6_i and hence P_i_. 

In fact it is easier to solve for g_i = —d^-i/dR. We do this by multiplying equation 
(46) by _R 3 (t -1 ), substituting equation (47) and differentiating with respect to R in order to 
obtain an equation for g_\ in the form 

R 2 g'U{R)+{?,--l)R g'^W+g^R) ( (3 - £)(3 - 2^)^^ i? 2 - 2(4 - 3^)) = 0,(48) 

where the prime indicates differentiation with respect to R. We have introduced the notation 
(r(x) denotes the usual factorial or 'gamma' function) 

I = r /2 ( CO s#) 2 (sin^^=(v^/4)r((g + l)/2)/r((g + 4)/2), (49) 
Jo 

p(j) is as defined previously, the polytropic index is n = p + 3/2, and we note that 
I 00 = 167ry(2)C 7 ^)+3/2) / 2p(f)+1 , (50) 

A QO = 16tt V / (2)C'7 ( ^ )+1/2) /2p(f)-i- (51) 

Here A 00 is the integral that appears in the integral over d 3 ( of dP_i/dR, namely 

A 00 = 8nV2C £ | eoo |(p-iy 7 _| eoo | rf| eoo |. (52) 

We have supposed that e DO < for bound particles. 

If we recall the definition (18) for the quantity 7 then equation (50) gives the explicit 
relation between I QO and C. 

The solution to equation (48) is (b 2 > 0) 

9-i = RWQV (d J m){2 - V {bR) + C 2 Y {3/2)(2 _ v (bR)) , (53) 
where 

if = (3-H)(3-2£)(^l)=2" (54) 
From this result and equation (46) we may find 0_i as 

4^0-! = -(1^) + ^^), (55) 

which should be compared to the zeroth order density from equation (17). 

We see by looking at the asymptotic forms of the Bessel functions near R = 0, that the 
second term in equation (53) varies as R~ 2 there, so that it allows for a point mass. We omit 
this by proceeding only with the first term, which in turn varies as near the origin. 

Figure (2) shows the variation of = O + «0_i for j = 1.125. We see that oscillatory 
fine structure is detected, which corresponds to the phase sheets seen in the early stages 
of the simulations (e.g. Merrall and Henriksen, 2003) before the the numerical smoothing 
instability. This shows that the analytic approach retains the phase sheet structure as indeed 
it must in the absence of any external perturbations. Of course to completely define the 
sheets would require higher order terms in the expansion. At this order the fine structure 
reflects the collective oscillation mode of relaxation. In the non-linear limit we expect these 
to become the distinct phase sheets. By contrast the coarse-grained result yields the full 
chaotic development due to the inevitable presence of numerical uncertainty. 



o n r\ a t~> A o tv tivtd A o r\r\r\ 1 0^7 



12 R.N. Henriksen 



0.001 



0.0008 - 



Theta 



0.0006 



0.0004 - 



0.0002 - 




Figure 2. Wc show an illustrative plot of the fine structure implied by equation (53) in terms of = © + o©_i as a function 

2 a 5 a 

of x. © is measured in terms of 7 00 /r T , x <— r/r oy and aC\r T /(AnGI 00 ) = 0.5 in order to show clearly the structure. The 
fiducial radius r is arbitrary. 



The limiting case f = 1 must be treated slightly differently. Equations (43) and (44) still 
apply with the appropriate and P , but we must remember that ( = Y = v. We proceed 
by calculating 

AnGQ 1 = AnG J P_ x d 3 v = A^^), ( 56 ) 

then taking the derivative with respect to R and using equation (43) for dP^i/dR to obtain 
the equation for g_i in the form 



R dRWdR {R j P 2 "" 



-9-i = 0, 



where now 



2nGI, 



\— [ v 2 exp-(v 2 /A-kGI 00 )AixC dv, 

jinn J 



and 

I 00 = 47rC J v 2 exp -{v 2 /AnGIoo) dv, 
so that 

A 00 = -1/(2ttG), 
and 

1 1 



(47r 2 G) 3 C 2 ' 

Equation (57) is now an Euler equation and the solution takes the form 



(57) 

(58) 

(59) 
(60) 
(61) 
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g-i oc R e , (62) 
where 

i = -1/2 ±iyjlji. (63) 
The general solution thus takes the oscillatory form 

g _ x = B J R- 1/2 sin( v /7/41n J R + 0), (64) 

where B and <fi are arbitrary. 

In this case the oscillations are less pronounced until the divergent centre is approached. 
This case is notoriously difficult to calculate in simulations so that it is interesting to establish 
the continuity with the polytropic behaviour found here. 



3 BREAKING THE SELF-SIMILARITY AS R -> 
3.1 Flattening the Singular Isotherm 

We break the self-similarity in our method by permitting an explicit dependence on t (Carter 
and Henriksen, 1991). We begin by showing the method in the algebraically simple case with 
j — 1. Returning to the section on coarse-graining we can write the zeroth order equation 
with time dependence as 

td t P - Rd R P - 2P = 0, (65) 
which has the general solution 

P. = (66) 

where u = Rt, and the same characteristics as for the first order term below. 
Proceeding to the first order term we find the equation 

-tdtP! + Rd R P 1 = -2P 1 + v R d R P - ^-—d v P , (67) 

dR v 

which has the characteristics 

dt dR „ . . 

* = ~ds ~ < 68) 

and hence 

^ + 2P 1 = g (-2P DO - V -d v P 00 + Rd R Poo) ■ (69) 
Here we define 

V ^ (70) 

for convenience. This quantity is AnGI 00 in the self-similar coarse-grained solution. Note that 
v and u (by direct calculation using the characteristics) are constant on the characteristics. 

If we measure t from the onset of the similarity breaking, then by continuity at t = we 
expect 

P 00 = C(u)exp-(v 2 /y) (71) 

since y — > 4^(7/00 = y(0) in this limit. We also require C(0) = C, the constant appearing in 
P 00 in the self-similar coarse-graining. The equation (69) becomes simply 
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j- g (P ie 2s ) = v R e- s Rd R P 00 = VR e- s P 00 (^- - v 2 Rd R (l/y)). (72) 

It is now necessary to write an equation for y(R, t) using Poisson's equation. From equation(16) 
we have easily 

±(Ry) = A^GC{u){Ryf"R^\ (73) 

which yields y = const = y(0), C = const = C(0) as one solution as is required at t — 0. In 
addition it has the general solution 

= _ 4 ^ G f 9m dR , _ Ki (74) 

where K is a constant of integration. 

Since we are interested in the origin of the flattening we proceed by holding C(u) = C. 
It is easy to see from equation (72) that a function C(u) that increases with u will assist 
any initial flattening, which is presumably due to the wave-particle relaxation itself. This is 
because a reduced central phase density (decreasing C with decreasing R) corresponds to a 
reduced space density for the same volume of velocity space. We are now able to write the 
solution for y as 

^l^^GC-KR 1 / 2 ) 2 . (75) 

We may note here that at a fixed r (with K > as we require below) 1/y increases with 
increasing time just as was found in the numerical simulations of Merrall and Henriksen 
(2003). This time-dependent dispersion is presumably due to the collective relaxation itself. 
Moreover at large r at any t, 1/y increases without limit which gives a cosmological DF 
peaked around zero velocity at large distances. Moreover in this same limit y oc 1/R, which 
implies a Keplerian potential outside most of the mass at large r just as was assumed and 
shown to give the r~ 3 exterior density profile in Henriksen and Widrow (1999). 
We now obtain from equations(72) and (75) the interesting result 

- {P ^)=e*v R v 2 P 00 [-^--y (76) 

We see that there will be a flattening if the first term in brackets dominates. This must be 
transitory in time since initially the two terms are equal (remembering equation (61), and 
finally y becomes large compared to the initial value. If is to be small then by equation 
(75) we must have 

K « ^17^' ^ 

where Rf is the radius where the flattening begins. We now approximate 1/y by the first 
term in equation(75) for R < Rf and integrate equation (76) using the dominant first term 
to obtain 

Pi » -(4ttGC)V^^. (78) 

Integrating this result in turn over velocity space (using cylindrical coordinates d 3 v = 
2nv ±dv±dv R ) we find 

01 = ~(Gcy (i6n3)w R °' (79) 

Consequently to this order 
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Figure 3. We show a plot of the scaled density, p/po against x = r/r and r = yJGp i according to equation (35) over the 
indicated ranges. We see that a core is beginning to form at the 'centre' of the system when r/(Ax) = 0.75 



on using equation (61). Now recalling that the physical density is e 2aT Q plus the definition 
of R we obtain the dimensionally correct result 



In this formula I 00 = p r 2 in terms of fiducial values in the self-similar region. It is clear that 
the flattening increases with dynamical age of the system at a given radius. 

We show in figure (3) an example of the type of flattening this equation implies. The 
figure shows p/p Q against x = r/r Q and r = \JGp Q t. We observe that by the time that t/Ax 
reaches 0.75 a central core is forming. 

There is no other physics in the problem but the collective relaxation provided by the 
wave-particle scattering. The time-scale of the density evolution seen in equation (81) is 
essentially a few crossing times at r. 

However we know from numerical work (e.g. Merrall and Henriksen, 2003) that suf- 
ficiently asymmetric initial systems will not be fully relaxed before the small scale time 
dependence is damped away, presumably by Landau damping. Thus the time that appears 
above is measured from each new epsisode of relaxation, that may be stimulated by mergers. 
For this reason and because of the demonstrated meta-stability of equilibrium polytropes, 
we believe that the relaxation can proceed through the polytropic sequence that we have 
found in (31). 

For this reason we discuss briefly in the next sub-section the similarity breaking for this 
sequence. 

We note that low mass systems are older in CDM cosmology and so they may be expected 
statistically to be more relaxed. This is the basis of a possible statistical test of equation 
(81), but this is best left to a separate paper. 




(80) 




(81) 



n n r\ a t~> A o 



TV IMT) A o r\r\r\ 



16 R.N. Henriksen 



3.2 Flattening the Singular Polytropes 

Since the method is precisely that of the previous section, we content ourselves here with 
giving a catalogue of the analogous formulae. We will discuss the results carefully however 
as they are somewhat surprising, due principally to the division at n = 5 (j = 5/4) between 
spatially finite and spatially infinite systems. 
The zeroth order scaled DF is found to be 

P = P 00 ((,u)R-V-f\ (82) 

where 

u = Rt^, C = yR {f ~ 1} , (83) 
and the coarse-grained characteristics are 

R = e- s , dt/ds = -t. (84) 
The equivalent of equation (72) for the first order scaled DF is 

+ (3 * - 1)P 1 = (~> - , (85, 

ds a i? 4 s y dmu e OQ J 

where by continuity at t = as above 

P o = C(u)\e 00 \ p , (86) 
and 

1 ^d* 

7 = 



R^f- 1 ^. (87) 



2(f - 1) dR 
Formally here 

e oo = y-70M, (88) 

just as in the self-similar coarse-graining, but now we must solve for 7 as in the previous 
section. For convergence we work as usual with bound systems for which e OQ < 0. 

Using the Poisson equation in zeroth order and proceeding as above we find that 7 
satisfies 

dw _ 32y/2V GC(u) {f) 

dR-^T hp+1 R^)l^) W ' ^ 

where w = "f(R, t)R^~ 2 ^ for brevity and n(j) is as in equation (33). 

When C(u) = C(0) = C, this equation has the solution w = ^R^ 2 ^ with 7 = 
(2ttGI 00 )/((j — 1)(3 — 2j)) as required. In general the solution may be written (assum- 
ing C constant) in the form 

where 

3 — 9- 

9=Vff- (91) 

5 l 

When j > 5/4 the system must have an 'edge' and so we set KR q E equal to the constant 
in the bracket of equation(90), where Re is the current edge of the system. Note that K > 
and that 7 — > constant as i? — > 0. 

When ^ < 5/4, the solution extends to infinity where 7 must go to zero so that we 
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must take K < in equation (90). We note that in this case there will be a radius given by 
KR q E ps bracket constant, outside of which 7 is dominated by the second term in the in the 
bracket. 

We are now able to perform the calculation indicated in equation (85) and so find 

Pi = f } Y f°r K^R^». (92) 

Integrating this over Y R cPY = R 4{1 ~^27i( R dC R C ± dC± yields finally the 'flattening' first order 
scaled density in the form 

Ox = ^ ^CK R ( 3+ *T^ } . (93) 

We remark that for the finite polytropes as R < Re 7 is constant, and since K > 
we have flattening as R — > provided that the power of i? in Oi is more negative that 
—2 j, which is the power of the zeroth order singular density. This requires j > (\/l3 — 
l)/2 ~ 1.303, which corresponds to n < 4.303. Between this value and | = 5/4 there is no 
asymptotic flattening as R — > 0, although there may be a 'shoulder' effect at R between R E 
and 0. Recalling that R = re~ ST = r / (at)~ , we see that such a shoulder region expands in 
time. 

When I < 5/4 there is a regime R > Re where 7 is dominated by the second term in 
(90). Then in equation (93) 7 is replaced by \K\ and the power of R becomes 

0! oc R- (3+p \ (94) 

and so there is always flattening in this regime, which is very rapid as j — > 1 since then 
p — > 00. This flattening will cease however as R — > Re from above, so that once again it 
is only detectable in this order in a shoulder region. This may be related to the fact that 
polytropes have a shoulder domain outside of the core where the local slope is less than — 2j. 
Thus we detect a flattening relative to the self-similar slope beyond this domain. The same 
effect is present in the isothermal case discussed above since 1/y ultimately becomes 'large' 
at small R and the flattening vanishes (1/y in equation (76) dominates the first term). 

We conclude from this analysis that the singular isothermal and polytropic spheres will 
tend towards the more general family of non-singular isothermal and polytropic spheres. 
This is supported both by the direct calculation of flattening found in this section as a 
result of the broken self-similarity, and by numerical simulations (Merrall and Henriksen, 
2003). Our analytic calculations conclude correctly that a steep region exterior to the core 
develops as part of this evolution. 



4 CONCLUSIONS 

In this paper we have focussed on the precise results that may be found by the resolution 
expansion method for spherically symmetric, isotropic collisionless systems. By insisting 
that the coarse-grained equilibrium structure should have no sub-structure (above the scale 
where the DF is well-defined) we show that the self-similar infall of collisionless matter 
assumes either a polytropic DF or a gaussian DF depending on its similarity class. These 
are compatible with the similarity density profiles that are strongly cusped at the centre of 
the system. 

However these singular polytropes (including the isotherm) are isolated mathematically 
in the family of mass distributions that are consistent with the polytropic DF. They are held 
in place so long as self-similarity is assumed by the pre-assigned 'class'. We have therefore 
studied the explicit breaking of the self-similarity that must occur at the centre of the 
system. This study shows that in lowest order the density profile tends to flatten relative 
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to the self-similar value, at least in a 'shoulder' region outside the true core. The flattening 
is progressive in time, and we conclude that this will be more evident in dynamically old 
systems such as dwarf galaxies, rather than young systems such as massive clusters. However 
the time-scale of the flattening is only several crossing times since it is due to the collective 
wave-particle relaxation, so that each system should be considered on its merits. Since there 
is a link between the statistical age of a system and its mass in CDM cosmology, this effect 
should be testable. We leave this to a less pedagogical paper. 

Our major suggestion is that this flattening, together with the evolution towards a central 
gaussian DF shown in numerical simulations (Merrall and Henriksen, 2003) for destabilized 
polytropes, should lead to the development of a central polytrope in collisionless matter. The 
polytropic index is given in equation (33) in terms of the similarity class. We expect this to 
be larger than 1 since both in principle and in practice the inverse case is unstable to rapid 
evolution towards 1. The index may be specified in terms of the initial perturbation that 
created the infall, or perhaps in terms of other dominant effects such as tidal truncation. 

This modifies slightly the theoretical expectations of Lynden-Bell and Nakamura (ibid), 
but this is probably because the polytropic state is only intermediary. Systems may stall in 
a polytropic state (or indeed one that contains evident sub-structure) since it seems that 
sufficiently asymmetric initial configurations cease relaxing collectively before a true gaussian 
is attained (Merrall and Henriksen, ibid). We believe that subsequent disturbances to the 
system such as those provided by mergers, can restart the collective relaxation and move 
the system farther along the road to equilibrium. This must be tested in detail however. 

We have pointed out both here and in Henriksen and Widrow (1999) that there is an 
inevitable edge effect in these systems when one is outside most of the mass. This 'law of 
the edge' that we also refer to as 'Keplerian' is simply p oc r~ 3 when the outer regions 
are dominated by a primordial density power law. However faster cutoffs are possible if for 
example tidal truncation plays a role (p oc r~ 4 in the limit), and indeed some such regime is 
essential if the halo is to be of finite mass. 

The fitting function that we recommend for the density profile of a not too evidently 
unrelaxed system (i.e. if the system were evidently bipolar one would consider each part 
separately) is shown in figure (1). The phenomenological parameters of other fitting functions 
are replaced here by the index of the central polytrope, since the fit to the law of the edge 
must take place where the slope of the polytrope is —3 (see figure). This index is in turn a 
direct link to the initial conditions of the dark matter infall. The degree of flattening is a 
direct link to the dynamical age of the system. 

The density profile for the metastable polytropes is found from the Lane-Emden equation, 
assuming that they flatten in accordance with their distribution functions inside the power 
law, self-similar range. The high order polytropes (n > 5) have to be terminated arbitrarily 
to avoid infinite mass, just as does the non-singular gaussian. 

It is interesting that in a recent paper by Roy and Perez (2004), wherein the evolution 
towards equilibrium of a collapsing system of N non-interacting particles is studied by means 
of a tree-code, they find that the equilibrium density profiles are fitted equally well by a 
high order polytrope (the Plummer sphere) or by a non-singular isothermal sphere. 

All of this behaviour seems to be a direct result of time dependent relaxation that is 
inevitable in such a system, just as was envisaged in the term 'violent relaxation' by Lynden- 
Bell (ibid) long ago. Of course the relaxation is effective without really changing the particle 
energies 'violently', and so we prefer to refer to it as 'collective relaxation' or 'wave-particle 
interaction'. This latter description evokes the Landau acceleration and damping that must 
be occurring. In any event we do not need physics beyond that of collisionless self-gravitating 
matter to understand at least the qualitative nature of the dark matter halo density profiles. 
The basic mechanism appears to be the retarded infall of the more massive shells for j < 3/2. 
The appendix shows that this effect is not restricted to spherical symmetry. 

In addition to these practical conclusions we have retained a number of pedagogical 
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Figure 4. We show on the left panel the density profile of an n = 2 polytrope, with an isothermal sphere superposed (always 
the curve with the larger asymptotic density) that has a velocity dispersion equal to 0.5 , f o . The densities are normalized to 
the central values. The right panel shows the same information for an n = 5 polytrope, with an isothermal sphere superposed 
that has a velocity dispersion equal to 0.2^ o . There is little distinction near the centre of the system. The scaling used is given 
in the text. The lower left panel shows an n = 9 polytrope with a superposed Gaussian having a dispersion equal to 0.1\t o 

aspects in the present paper. We retain the mathematical details in order to be as convinc- 
ing about the utility of the resolution expansion method as possible. In addition we have 
included a section on fine-graining, where we show that in lowest order the phase space 
structure found in numerical studies is partially recovered. This sub-structure appears as 
outward propagating waves in linear order (R = r/(ai«)), which we take as an interesting 
confirmation of the 'wave-particle' nature of the interaction. In any case it shows that the 
expansion method is rather flexible and reliable. It is in fact more general than in this appli- 
cation to self-similar infall, since any parametric non-canonical transformation may be used 
(e.g. paper I) 

We show in figure (4) for reference the density profiles for the n = 2 polytrope (corre- 
sponding to a constant angular momentum including zero and j = 2) and for an n = 5 
polytrope (a Plummer sphere). Isothermal spheres that fit well with the polytropes at the 
centre of the system are superposed. We have measured the densities in terms of a common 
central density p Q and we have used the same scale factor b 2 = l/(4irG^- 1] c n ), in the 
notation of Binney and Tremaine (ibid). This requires \l/ = 47rGp b 2 and = p /c n . Thus 
with p and b chosen the curves are determined physically. The remaining parameter in the 
equations is ^ /o~ 2 , and we show the isothermal curve for ^ /a 2 = 2 together with the n = 2 
polytrope and that with \l/ / a 2 = 5 together with the Plummer sphere. The profiles are not 
easily distinguished near the centre of the system. When n = 9 the curves can be made to 
fit over a wide range as shown. 

In an appendix we present the ambiguities that arise in a study of a completely general 
system. There is always an access path to the polytropic/gaussian state however, and general 
arguments due to Lynden-Bell, Nakamura (ibid) and Aly (1989) allow similar conclusions 
to those found for this special case. 
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Appendix: Collisionless Systems in 3D with no imposed symmetries 



In this appendix we contrast the behaviour of an accreting halo that is permitted to fill 
a full six-dimensional phase space, with the restricted halo studied in the body of the paper. 

The relevant equations are those of section (2.1), before the restriction to the spherically- 
symmetric isotropic case. 

We proceed to the expansion in terms of reciprocal powers of a as usual, but we note 
from the variable transformations that this represents a coarse graining of phase space only 
if j < 2. Should | = 2 we have an interesting special case where there is a global integral 
having the dimensions of angular momentum. Moreover the transformation is canonical in 
this limit so that there is no change in resolution as a increases. We may still expand in a 
series in a" 1 however and this corresponds (see the time transformation in5) at a fixed T 
to a cruder resolution in time (At = e aT AT). We know from Fillmore and Goldreich (1984) 
and Bertschinger (1985) that if the halo develops from a primordial density perturbation 
characterized by an effective power law oc r~ e , then the similarity class is given by (36) in 
the text. 

Since we expect e to be positive there is no value that gives f = 2. Thus we must abandon 
setting the initial density profile independently of a limiting specific angular momentum J. 
With Newton's constant G this determines a density profile (38). This may well be relevant 
to halos formed from material that has been spun up tidally (presumably that which has 
fallen in from larger spatial scales). The outer regions of globular clusters and of elliptical 
galaxies may reflect this behaviour in their density profiles. 

That specific angular momentum is not always conserved in spherically symmetric po- 
tentials (see below) is due to the time dependence of the potential in these non-isolated 
systems. The mass growth inside a given elliptical orbit is not adiabatic. 

Proceeding with the expansion the zeroth order DF may be found by the method of 
characteristics to be (we drop once again the explicit time dependence, which is the condition 
for self-similarity) 

Po = Poo((,0A,Ro)e- {3 «- l)s , R = R e- S , ( = YR^- l \ 



e, 



/ Pod% (1) 



where ( and 9, 0, R Q are constants on the characteristics. We suppose for simplicity that 
all characteristics cross the same R = R Q so we may set this equal to 1. This is most likely 
to be a (scaled) radius that each particle crossed on entering the system since it is clear 
subsequently that particles are trapped near the centre of the system and we are forbidden 
to take R Q to be arbitrarily small. It is clear that this choice determines the origin on each 
characteristic although each particle is in general present at s = at very different times. 
In this case we have 

s = jlnR, (2) 
and the zeroth order DF as 

P o = Poo{l0^)R { ^ ] - (3) 

The formula for the scaled density becomes 

e = ir 2 ? / p 00 d 3 c = ioo(o, 0)/r 2 l (4) 



It is worth emphasizing here that whenever the series is stopped at the zeroth order we 
have a steady solution; since by definition ( = ■ur < -t~ 1 ) ) which is time independent along 
with 9 and 0, and since from the various definitions 
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fo = P o e^ T = P oo (C,0,<P)r- {3 ^ ) . (5) 

Another way of looking at the result is to note that a is effectively infinite when the series 
is exact at zeroth order, and thus there can be no time dependence since such an a implies 
an average over all time. 

To continue the expansion to higher orders we must consider the solution of the Poisson 
equation in zeroth order. This is 

V 2 m = AnGR-^I 00 {d^). (6) 
One form of the solution in spherical harmonics is 

*.= ^i2< 2 - 2 *U(M), (7) 

l \8> 

where = (2 — 2|)(3 — 2|). Here the function A is given by 
A->: A lm Y lm (0.&). A lm - r ~ ! ' 



i(f)-/(/+l)' 

iL = I hoYL dn. (8) 

The spherical harmonics Y\ m are used with the convention adopted in Binney and Tremaine 
(1987) and must not be confused with the scaled velocity used above. The superscript * 
indicates complex conjugate and £, m are the usual integers. The value £ = 2 — 2j must be 
excluded, and we note for future reference that A/i(j) is finite and negative (£ = 0, £ = — 1 
are excluded) asi-^0or|-^3/2. 

The procedure is slightly different in the excluded case of j — 1. A sum separation yields 
a general result (a product separation requires a proportionality between A and I 00 ) in the 
form 

^ = k\nR + A{9.<j>), (9) 
where A(9 ) 4>) satisfies 

-±-d e {sin6d e A) + -^dU = AnGI 00 - k. (10) 
sin V sin V v 

Here A; is a true constant. 

Proceeding now to the first order in the expansion we obtain for the next term in the 

DF (taking care to differentiate before evaluating on the characteristics (such that ( is R 
dependent) we find that 

ds v ' 

R~ 3 {-(3 - j)CrPoo + l)C R Cd,Poo + (ed e P 00 + ^-d^ + ((e + Cj)^P 00 

s Sill u 

-2(1 - -)(^GA/i(-))d (R P 00 + (CJ cot^ - CnCe - ±nGd e A/i(-))d (9 P 00 

-((4,(r + Uocot(9) - i{ a^ n{d f ^)dc,Poo} ■ (11) 

This result is of interest in itself since, as has been discussed in paper I, it can be used to 
describe the deviation from the self-similar density law near R = 0. However for our present 
purposes setting this expression equal to zero yields a linear partial differential equation 
for the exact steady P OQ in the intermediate self-similar domain as in the examples above, 
and this may in turn be converted to the following holonomic system of equations using the 
method of characteristics: 

sts onr»/i d A o attvto a c r\r\r\ 1 o 1 ? 
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d( R a 



1)CI + Ce + Cj" ^GAj (3-2- 



^ = (|-l)C^ + C|cot(^)-C^-47rG'^/i(|), 

^ = (|-l)C^-C^cot(^)-C^-47rG^(^)(sin^), (12) 

de 

da 

d<j> 

da 



= Ce, 

= C*/sin0. 



For the excluded case with f = 1 the same procedure using the zero order potential of 
equation (9) yields (we write them again for convenience but note that it suffices to replace 
4nG/i by 1 and AnxGAj (3 — 2^) by k in addition to j — 1 to obtain the following) 



= Cjcot(0) - CnCe - d e A, 

= -C*Cfl-C*C*cot(0)-fyMan0), (13) 

= Ce, 

= C^/sinfl. 



dPoo 
da 

d ^ R c 2 4- C 2 - k 
da ~ *' 

dQ 
da 
d(<t> 
da 
d6 

da 

d(f> 

da 

Here a is taken to be the 'path parameter' in the phase space of this equation, but it will 
normally not appear in physical quantities. 

These equations may be rearranged in two important ways. In the first instance we can 
combine the third and fourth of these equations to obtain 

I<I = _( 2 _ ®)( R (l - AnGdA (14) 
2 da 5 i da 

whence on eliminating using the first of these equations one finds 

rfln(P 00 (Cj)^) 4vrG(3 - f ) \_d_A 

da ~ <(f)(2-f) CI do' 

where we have used the notations 

Ci-Cl + Cj, <) = ^fy- (16) 

One simply sets A^G/i = 1 in these equations in addition to j — 1 in order to obtain the 
result for the excluded (isothermal) case. We observe that there is never an integral unless 
dA/da = 0, which holds strictly only in spherical symmetry . 

However and more importantly, considering for the moment the general case, the last 
five equations of (12) may all be combined to give 

d In f^lli + ^ga\ 

U 111 ^ „ T ;/a\ ; q, 

= ~ 2(1 - J** (17) 



(15) 
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and this in turn may be combined with the first equation to write the the steady, stable, 
self-similar DF in the form 

P 00 = S(« i )|e|^. (18) 

Here Ki denotes any other integrals, the index p(j) is as in equation (2.2)) of the text, 
and the scaled specific energy e is 

(Ci + C? + Cf) , ^GAje^) 

2 + i(f) • {W) 

If we now reconstruct the physical DF from the various definitions we find that 

/ = a(«o \e\«v, (20) 

where the specific energy is 

E= {VI + V 1 + ^ + ^r A <Pl (21) 
if the potential is 

$ = ^GA{6A) r{2 - n) _ (22) 
Hf ) 

In general any will be constructed from quantities invariant on the characteristics of the 
CBE so that Jeans' theorem is satisfied. We have not in fact found any other explicit integrals 
in the general case as might be expected in the absence of any phase space symmetry. If 
such integrals exist they are time independent since all of the variables in equation(12) are 
themselves time independent. 

We should note that this DF is the only form that is steady in zeroth order independently 
of the higher orders and it may in fact be deduced this way (cf HW99). However the present 
procedure shows explicitly that it is a stable self-similar form that may be compared to the 
isotropic result. It is more general in terms of the gross asymmetry that is permitted. In this 
case the form is not obviously unique, assuming that other integrals exist. In the absence of 
any imposed symmetry however it seems likely to be relevant with B = b a constant. This 
constant has to be determined self-consistently with the Poisson equation of course. 

Proceeding in the same fashion with the last five of equations (13) for the case j = 1 we 
find directly that 



2 e 



P 00 = B( Ki ) e~— , (23) 
where the scaled specific 'energy' is (the radial part of the potential does not yet appear) 

e=&±£+A(6,<l>). (24) 

This leads through the various definitions to the steady, stable, self-similar DF in this 
isothermal case as 

f = B(Ki)e-™, (25) 
where the specific energy is 

E= V llA + + A(9, ( p) + k\nr. (26) 

Once again we do not find any other explicit integrals, and since this is the gaussian form 
expected on theoretical grounds (Lynden-Bell 1967, Nakamura, 2000) we suspect it is the 
unique form. It can also be deduced by requiring a steady DF for j — 1 without the present 
formalism. It is an important result in any case as it shows that the density profile (1/r 2 ) 
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that appears so frequently in numerical simulations is also associated with the gaussian DF 
expected on theoretical grounds (ibid) in general symmetry. 

It should also be noted that we have omitted discussing the final phase of any specific 
choice of DF, which is to solve the resulting Poisson equation for the potential function 
A(0,<f>). If we consider only what we suppose to be the limiting case when j — 1, then we 
have to solve the isothermal type potential equation (10) with 

ho = e — j B{Ki) e~~ d 3 (. (27) 

Even with B constant, we see that there may be only certain values of k that are per- 
mitted (for example k = lim^o An G loo in order to avoid a singularity on the axis). One 
solution is clearly A = constant and k = 4ttGI 00 , that is spherical symmetry, which is 
clearly a minimum energy condition for positive A. 

The special case j — 2 is also of interest as remarked above. The equations (12) now 
yield 

f = B(K 1 K i )\E\ l ' 2 ) (28) 
in agreement with equation (20), but now in addition we have the integral 

C 2 

K = + 2nG A(0,<j>), (29) 

which certainly removes the strict uniqueness of the DF. The variables ( here have the 
dimensions of 'actions', namely oc rv. The radial dependence of the potential is r~ 2 so that 
we are strongly forbidden to continue this case toward the centre of the system. Nevertheless 
this case illustrates nicely the fact that with more 'isolating ' integrals the DF is less unique. 
It seems that it is best defined in the absence of integrals in addition to the energy, when all 
of the particle orbits would be 'irregular' and would tend to fill the available phase space. 
This would mean that the Jeans theorem is irrelevant if not wrong (cf Binney, 1982). 

We proceed to discuss the special cases of spherical symmetry in physical space (not 
necessarily isotropic). Such limits are important since we expect for example the spherical 
isothermal case to provide the accretion equilibrium for all 'flat' cases (f < 1) and since 
spherical spatial symmetry in general has been claimed to minimize the energy of collisionless 
systems (Aly, 1989) subject to a reasonable phase space volume constraint. 



1 Spherical symmetry with anisotropic velocities 

With the assumption of spherical symmetry in real space but not in velocity space, A is a 
constant and equation (15) yields immediately; 

P 00 = B(K l )(tl)-^\ (30) 

where the Ki denote any other integrals. In fact the third, fourth and fifth of equations (12) 
may be combined to yield another integral in the form 



, =((!) BfaM + ^. (31) 



6- 



This appears to be the only other explicit integral. Using the various definitions we obtain 

k = (j!)^E, (32) 
where j\ = (rv±) 2 and the specific energy is 

E , *1±A + !!E^V- 2 f>. (33) 

onn/1 D A O A yTTVTTD A C C\C\C\ 1 O 1 ? 
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Consequently we find the steady self-similar DF in the form 

f = B{K,K l ){jl)-^\ (34) 

in agreement with Jeans' theorem once again. But we know that in general the DF has the 
form of equation (20) and we find this to be so if we take B(k, «j) to be 

S(«,« i )= J B(« i )x(|«|) 5 ^, (35) 
for then the DF above takes the general form 

(3— a ) 

f = B(Ki) x \E\**% , (36) 

where now k itself may be one of the Kj. Thus the DF of equation (36) is clearly not unique 
in spherical symmetry. In the case of j — 2, the form is given correctly by equation (36), but 
the additional integral is The inhibiting effect of angular momentum on the relaxation 
was already emphasized by Aly (1989). The DF (36) was already found in Henriksen and 
Widrow (1995) as part of a family of strictly steady self-similar solutions. 

We remark here that the preceding results are those that would have been found in paper 
I for this case, had an algebraic error not occurred there. It is easy to demonstrate this fact 
by applying our procedure directly to the spherically symmetric form of the CBE as in paper 
I namely 

d t f + V r d r f + - d r $)d v J = 0, (37) 

where the tangential angular momentum is jj_ = r 2 v\ and all derivatives are taken while 
holding it constant. The procedure followed above leads directly to the conclusion of equation 
(34). This corrects the result in paper I and confirms the treatment given here. The result 
with j — 1 is given in the next section. 

In spherical symmetry in an isolated steady state we know that the three components of 
angular momentum are also isolating integrals, and we might expect these to appear in the 
ultimate DF, even though only jj_ appears explicitly in our analysis. This is because of the 
non-isolated nature of the system as it grows which leads to the form (36). We note that 
regularity in jj_ suggests that k does not appear again in equation (36), if the function B is 
restricted to be a power law. 

Thus our major conclusion in this section is that with only spatial spherical symmetry 
we can not determine uniquely the DF by looking for stationarity and exact self-similarity 
alone, unlike the result for the one dimensional systems studied above. Should there be no 
additional 'strange' integrals, and if we require regularity in j± together with at most a 
scale- free function B in equation (34), then we re-acquire a polytropic form (equation(36) 
with B constant). 



2 Spherical symmetry: f = 1 

This example is important because it is expected to be the state towards a system ultimately 
tends while relaxing, and it is not necessarily isotropic. We find on putting A = equation 
(15) that 

P OQ = B(k) CI 2 , (38) 
where again we find only one additional integral (£ = v here) 

n= V -l±^- k -lnvi, (39) 
and so from the various definitions above with j = 1 we obtain the equilibrium DF as 
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f = B(K,K t )(rv ± y 2 . (40) 

We remark that the only choice for B(K,Ki) that removes a logarithmic singularity in 
the density (recall that d 3 v may be written ndvj_dv r in spherical symmetry) both as 

v 2 ± -> oo (41) 

and as 

v\ -> (42) 
is 

B(K,Ki) = B{m) x e~T?. (43) 
Hence we obtain finally from equation (40) the well-behaved DF in the Gaussian form 

/ = %)xe"T (44) 
where the specific energy is 

v 2 

E=— + k\nr. (45) 

This result now agrees with the general form as found in equation (25) but again it is not 
obviously unique since one of the /% may be k. Since however k is not well defined at the 
velocity limits, and since we have found no other explicit integrals, we suspect this form to 
be unique. It is also in agreement with the argument from maximum entropy (Nakamura, 
2000). 

Assuming this uniqueness, we may replace B in equation (44) by a positive constant b 
and so completely determine this result in physical terms. Thus we re-introduce a fiducial 
density and a fiducial radius p Q , r Q , so that from the Poisson equation we have 4nGp = k/r 2 
and hence 

k = 4nGp r 2 . (46) 

Then from the integral of the DF in (44) over velocity space (making the usual extention to 
infinite velocities) we have 

K> 3/2 >W7r (47) 

Consequently we conclude that in spherical symmetry, the isothermal (j = 1) steady 
self-similar DF is likely to be a gaussian with the above properties. This is much as in the 
isotropic case examined in the text. Since theoretical arguments referenced above reveal this 
DF as a maximum entropy condition, and since we might expect spherical symmetry to be 
also a condition for equilibrium, we suggest this form as a limit towards which collisionless 
accreting systems tend in many cases. We can not conclude this to be universally so because 
we have seen that steep initial density profiles (actually f > 1 which is more general) are 
prevented from reaching this state, apparently because of 'trapped' surfaces in phase space. 
However we note also that it is difficult to distinguish the density profile of a large index 
polytrope from that of the isothermal sphere as is demonstrated in figure (4). Our suggestion 
is that an actual system may be found in a polytropic state with an index (such that f > 1) 
that increases in time under the continuing excitation due to mergers until it is virtually 
indistinguishable from a (lowered) gaussian. 



